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Abstract 

Several viscous incompressible flows with strong 
pressure interaction and/or- axial flow reversal are 
considered with an adaptive multigrid domain 
decomposition procedure. Specific examples include the 
triple deck structure surrounding the trailing edge of a flat 
plate, the flow recirculation in a trough geometry, and the 
flow in a rearward facing step channel. For the latter case, 
there are multiple recirculation zones, of different 
character, for laminar and turbulent flow conditions. A 
pressure-based form of flux-vector splitting is applied to the 
Navier-Stokes equations, which are represented by an 
implicit lowest-order reduced Navier-Stokes (RNS) system 
and a purely diffusive, higher-order, deferred-corrector. A 
trapezoidal or box-like form of discretization insures that 
all ma« conservation properties are satisfied at interfacial 
and outflow boundaries, even for this primitive-variable 
non-staggered grid computation. 

Introduction 

Viscous interactions are typically associated with 
turbulent or high Reynolds number (Re) laminar flows. 
These interactions are quite frequently characterized by the 
appearance of high flow gradients that are most significant 
in small or ’ thin ' domains of finite extent, and in one or 
more directions, e.g, boundary or vortical layers/regions, 
triple deck structures, shock wave structure. Outside of 
these regions, the flow is generally more highly diffused or 
inviscid so that the flow gradients are less severe. 
However, the flow character in these smoother regions, 
which generally encompass a major portion of the flow 
domain, can be significandy influenced by the interaction 
with the high gradient layers. In order to accurately assess 
this of viscous interacting flows with discrete 

computational methods, (1) local grid refinement is 
required in the high shear layers, and (2) simple, efficient, 
adaptive methods, that effectively communicate information 
between the disparate flow domains, and at the same time 
maintain all conservation properties, are necessary. 

In the present investigation, an adaptive, multigrid, 
domain decomposition strategy is combined with a 
pressure-based form of flux vector discretization in order to 
a cco m p lish these goals 1- *. The governing Navier-Stokes 
equations are evaluated through an implicit, lowest-order 
in Re, reduced Navier-Stokes (RNS) subsystem 14 , that is 
combined, when necessary, with an explicit purely diffusive 
deferred-corrector (DC) in viscous layers. Local directional 
refinement that is driven by specified flow parameters and 
accuracy limits is achieved by sequentially splitting the 
overall flow domain into a variety of subdo mains . In the 
present analysis, this domain decomposition strategy is 
applied, in conjunction with an adaptive multigrid 


algorithm, in order to achieve the appropriate level of grid 
refinement. In this approach, each grid in the multigrid 
hierarchy, is of equal or lesser extent than all of the coarser 
predecessors. The subgrids are split into several 
multidimensional subdomains that are defined by specified 
directional and global resolution requirements. A similar 
approach has been presented for cavity and backstep 
geometries in a recent publication 5 ; although, no attempt 
was made to meet the differing needs for refinement in two 
or more coordinate directions. In the present 
investigation**, this is achieved with a subdomain procedure 
that allows for segmentally varying grid resolution in two or 
more directions throughout the flow field. This leads to 
more optimal grid refinement, and, through the adaptive 
multigrid procedure, information is very effectively 
transferred between high and low gradient domains that 
have distinctly different grid structure. In addition, the 
equation solver can differ from subdomain to subdomain, 
e.g., direct solvers can be used in strong interaction 
domains, line relaxation in moderate interaction domains, 
etc. 

In the present analysis, several two-dimensional, 
steady, incompressible, large Re laminar and turbulent flow 
examples are reviewed and the results are compared with 
other computations or experiments. The problems to be 
discussed include, the laminar trailing edge (triple deck) 
flow past a finite flat plate, the laminar recirculating flow 
associated with a trough geometry, and the laminar and 
turbulent flows in a backstep channel. The use of pressure- 
based flux splitting and a trapezoidal or box-like 
discretization for the implicit RNS subsystem leads to a 
precise prescription of the surface normal boundary 
conditions on the local subdomain boundaries. This 
ensures that interfacial and global mass conservation 
requirements are automatically satisfied. This is generally 
not the case with some characteristic-based Navier-Stokes 
schemes, where special conditions are required in order to 
satisfy interfacial and global mass conservation. . The 
primitive variable system considered herein is also directly 
applicable on non-staggered grids. This differs from many 
other incompressible primitive variable Navier-Stokes 
formulations, that require pressure Poisson solver or 
artificial compressibility concepts. 

(Inventing Equations an d Discretization 

The governing Navier-Stokes equations, shown here in 
sheared cartesian coordinates, is written for incompressible 
flow in non-conservation form: 

U{ + v, - 0 Continuity, (la) 

uu { + uy b (V+y b u) { + V[ (l+y b 2 )u,+y b V b l 

+ p? - _Lu^ + dc 5-momentum, (lb) 

u(V+y b u) t + V(V+y b u), + p, - DC n -momentum (lc) 


f . 

where £ = x ; n = y - y b <*> ; V - V - y^u is the 
contravariant velocity component in the q or normal 
direction (for y^x) < < l) ;y b (x) is the surface definition 
and (u,v) are the cartesian velocities in the (x,y) directions. 


difference for incompressible flow, so that the elliptic 
acoustic interaction or upstream influence is introduced 
through the p { * (pj+i - Pi)/(A£)i contribution. For 
reverse flow regions, additional negative eigenvalues or 
upstream influences appear through the convective terms. 


For turbulent computations, the k-e model is 
employed. This introduces two additional equations for k 
and e. These equations in cartesian coordinates and non- 
dimensional form are given as: 
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where v t - and c, f c c1 , c, 2 , a, and ff k are 

predetermined dimensionless constants which have the 
values 0.09, 1.44, 1.92, 1.0 and O respectively. A modified 
k-e model was recommended by Thangam [17] in which the 
value of the constant c^ is changed from 1.92 to 11/6. 
This was shown to give much better results for the backstep 
channel geometry. 

A three layer law of the wall is used at the upper and 
lower walls. This is given as: 

y* tor y'i 5 

u* * ■ -3.05 + 51ny* for 5 < y* £ 30 • 

5.5 + 2.51ny* for y* > 30 

where y* = yu,/v and u* = u/u,. 


A simple line relaxation procedure is used to solve 
the system of equations for k and e. The differencing used 
for the k,e equation is consistent with the pressure flux-split 
discretization. The law of the wall provides the boundary 
conditions one point away from either wall, which can then 
be used to implicitly to solve for k and e at any given 
station. The k-e equations are decoupled from the 
governing RNS equations. Each pressure relaxation sweep 
is followed by a sweep to solve the^k,e equations using the 
latest available values for u, v. This procedure is 
convergent. 

Grid Structure 

In general, the N" 1 multigrid level consists of several 
subdomains. Each multigrid level has an equal or lesser 
extent than the coarser grids of the multigrid hierarchy. 
The first two grid levels cover the entire computational 
domain. The mesh size is initially quite coarse in the 
directions in which adaptivity is to be prescribed. Each of 
the multigrid levels comprise several subdomains, which 
derive part of their topology from the subdomaini ng p attern 
of the coarser predecessor. Within each subdomain, of a 
given multigrid level the refinement is specified 
independently. Thus, each subdomain of a multigrid level 
can act as a parent for a subdomain or subdomains at the 
next finer multigrid level. If at a given multigrid level, a 
particular subdomain is refined in only one directi on, e.g., 
q, then on subsequent multigrid levels, further refinement 
within this subdomain is performed only in the q-direction. 
A similar strategy is adopted for the {-direction. Only 
subdomains that result from refinement of a parent 
subdomain in both the £ and q directions require further 
decomposition according to the direction selective 
refinement specifications. 


The RNS approximation is given by thelpwest-order 
system obtained by omitting the purely diffusive deferred- 
corrector (DC) terms. These terms are retained selectively 
in some subdomains, when they are important. The RNS 
system is in effect a composite of the Euler and 2nd order 
boundary layer equations 1 * 3 . Trapezoidal or ’box' two point 
(ij± 1/2) differencing is used for all normal derivatives in 
the first-order (in q) RNS equations (la,lc), and three 
point central (ij) differencing (in q) is applied for the axial 
momentum equation (lb). If a shear factor o = pUy is 
introduced in (lb), a box-scheme can in fact be developed 
for the entire RNS system (1). All axial (£) convective and 
pressure derivatives are upwind differenced with a 
pressure-based form of flux vector splitting 3 , wherein the p t 
term is represented, for compressible flow, by 

P* * «|-l/j(Pl"Pt-l)/^« + ( 1— <■*! *1/2) (Pi»i~Pt) /^< f 
and ■ [ym|/(1 ♦ (Y -1 )***) • 

Here, M. is the streamwise Mach number and y is the 
ratio of qjedfic heats. This reduces to a simple ’forward’ 


_ Refinement Strategy 

In most adaptive gridding methods, on any grid level, 
an estimate of the truncation error of the discretized system 
of equations is used to identify those regions that require 
finer grid resolution 3 . The overall truncation error 
estimates, however, do not provide information on the 
specific directions) that require refinement. Therefore for 
regions requiring higher resolution, the grid is refined in 
both directions, even though only one coordinate gradient 
may be significant In order to achieve directional 
refinement adaptivity it is necessary to monitor the 
truncation error of selected gradients or derivatives. For 
the problems considered herein, the truncation error for the 
pressure and vortidty gradients, e.g., p { and u„, are 
monitored in order to define the regions that require 
refinement in, { and q, respectively. Additional gradient 
parameters can be added when necessary. 

The truncation error estimate is obtained from the 
solution on two successive grids of the multigrid hierachy. 
In order to determine the truncation error in a £ (and/or 



n ) derivative, the finer of the two grids must have regions 
that are refined in the £ (and/or q) direction(s). Although 
the p t and u„„ terms are the key derivatives for the 
present analysis, the truncation error of these terms alone 
will not suffice to ensure that uniform accuracy is achieved 
throughout the flow domain. The global truncation error 
for the full discrete system of equations is monitored for 
this purpose. 

Two types of adaptive calculations are performed for 
the geometries considered herein. 

a. One-dimensional adaptive calculation (semi- 
coarsening multigrid), with adaptivity in the 5 
direction and with a preset stretched q grid. 

b. Two-dimensional adaptive calculation, in which the 
refinement is automated in both directions and 
unifo rm grids are used in each subdomain. Grid 
stretching is not applied, except as the grids change 
discretely from subdomain to subdomain. 

The underlying procedure is identical for both 
methods. The solution is first obtained on a coarse grid, 
for those direction(s) in which adaptive multigrid 
refinement is to be considered. For the semi-coarsening 
adaptive calculation, refinement is performed only in the £ 
direction. The grid is refined over the entire domain, and 
an improved solution is obtained. From the two full grid 
solutions, the truncation error of the key derivatives and 
also of the global discrete system is estimated using 
Richardson extrapolation. Two types of refinement criteria 
are used. In one procedure, a tolerance is set for the raw 
truncation error and, in the other, a tolerance is set for a 
truncation error normalized with the maximum value. The 
results obtained with the two methods, i.e., identification of 
the regions that require refinement in the respective 
directions) are quite similar. 

For the one-dimensional (in 5) adaptive calculation 
only one subdomain results. This decreases in extent as the 
grid level increases. For the problems considered herein, 
the significant flow gradients in 5 are centered around the 
small region [ £ | < (g. For more complicated flows, it is 
possible that disjoint subdomains in £ will result. For the 
two dimensional adaptive calculation, however, the various 
regions will have different refinement requirements; 
therefore, it is necessary to define regions that have 
disparate grid requirements. Subdomains requiring 
refinement in the q direction, or the £ direction, or in both 
(£,q) directions, are identified. Although different meshes 
are used in different regions, within each subdomain, a 
uniform grid is specified. This procedure is applied on the 
thir d and higher levels of the multigrid hierarchy. The 
calculation proceeds with intergrid multigrid transfers. On 
convergence, the truncation error estimation process is 
repeated with the N*^ multignd and the stored (N-l) 
multigrid level grid solutions. 

Multi grid Implementation 

For the RNS system of equations (1), without DC, a 
semicoarsening multigrid procedure has been presented 
previously 9,10 to accelerate the convergence of the global 
pressure relaxation procedure 1,1 . A von Neumann analysis 


of the linearized form of the RNS system shows that the 
rate of convergence of the global procedure is dictated by 
the marirmiTn eigenvalue X, as given by 

X " 1 - c^fAO^/q*, 

where c, is a constant of 0(1); N { is the number of stations 
in the £ direction; q M is ^e normal boundary location, and 
A £ is the axial step size. The convergence rate is 
significantly improved if the extent of the domain in the 
two directions is reduced. The ament multigrid domain 
decomposition procedure, in effect, reduces rj M whenever a 
fine A (is specified, and thereby achieves comparable 'coarse 
grid' convergence rates on fine grids. 

In the present application the multigrid method is 
implemented in a Full Approximation Storage (FAS) mode. 
The global pressure relaxation procedure considered herein 
essentially reduces to a block SOR procedure (in £) for the 
pressure in attached flows and for the pressure and 
velocities in reversed flow regions. At each station, an 
implicit, fully coupled tridiagonal system is inverted. When 
highly stretched grids are used in q to resolve the boundary 
layer, the semi-coarsening mode of the multigrid method 
has been shown to be more effective than the standard full 
coarsening mode. In this mode, the streamwise grid alone 
is coarsened when the calculation shifts to coarser grids. 
The same q grid is retained throughout Significant gains 
in the overall effort have been achieved with this 
approach 9,10 . 

A source term (1ST), first introduced by Israeli , is 
required in order to achieve satisfactory performance of the 
multigrid procedure 10 . The 1ST acts as a form of under- 
relaxation or smoother for the pressure field. This leads to 
much smoother residual fields, which are essential for good 
representation on the coarser grids. However the 1ST leads 
to a slower asymptotic convergence rate on any given grid. 
The do p™'" decomposition procedure reduces this 
limitation to a large extent Since the truncation error in 
the p, term is used to determine regions needing 
refinement in the £ direction, subdomains in which the grid 
is only refined in the q direction, will generally have a 
reasonably converged pressure field from the coarser grid. 
Thus it is possible to perform the multigrid calculation 
without the 1ST smoother in these subdomains. 

In the present investigation, the one-dimensional 
adaptive calculation adds an element of sub-domaining to 
the semi-coarsening analysis 9 ' 10 , so that only portions of the 
global domain require fine grid resolution in the £ 
direction. For the two-dimensional adaptive calculation, 
the multigrid algorithm is implemented in the standard full 
coarsening for domains that require refinement in both 
directions and the semi-coarsening mode for those domains 
requiring refinement in only one direction. One fine grid 
work-unit is comprised of one sweep in each subdomain 
belonging to a given multigrid level. This also includes the 
interdomain transfer processes. The decision to move the 
calculation back to a coarser grid is based on the rate of 
convergence on each subdomain. If the ratio of the 
residual norm between two successive global iterations, in 
any subdomain belonging to that multigrid level, falls below 
a certain value, typically 0.85-0.95, then the calculation is 
restricted to the coarser level. The fine grid solution is not 



corrected until the residuals in the coarse grid subdomains 
are all reduced to a value one order-of-magnitude lower 
than the ma ximum residual over all subdomains in the finer 
level. The multigrid components are summarized as 
follows, 

a. Relaxation: - s k u k _,, where S k is the global 

pressure relaxation operator and u k on convergence 
satisfies L k u k = f*. Here k represents the present or 
finest multigrid level and n represents the iterate. 
LV - f* is the discrete approximation of the 
continuous problem Lu = f 

b. Restriction to coarse grid where the following 

equations are solved: L k * 1 u k * 1 « l£ ’r k + L k ‘'i£ u n 
for points on the coarse grid which lie within the fine 
grid and L k " 1 u k ' 1 = f k '’ for points on the coarse grid 
that lie outside the extent of the fine grid. Here 
r k , f* - L k u k . ij' 1 and l k "’ are fine to coarse 
grid transfer operators. The full-weighting operator 
recommended by Brandt 12 is used to transfer the 
residuals and the solution was restricted by using a 
simple injection operator. 

c. Prolongation or correction where the fine grid 
solution is corrected with the solution from the coarse 
grid modified problem. 

' where x £-i * a 
coarse to fine interpolation operator 

It should be noted that in the present calculation, the 
multigrid tr ans fer operations play a dual role. In addition 
to accelerating the convergence of the relaxation procedure, 
they also tr ansmi t information from the finer grids to the 
coarser grids, and thus improve the accuracy of the solution 
in regions of the coarser grids where refinement was not 
required. The second term in the multigrid restriction 
process, acts as a truncation error injection term and 
improves the discrete approximation on the coarse grid. 
Thus on the coarser grids, instead of solving L ’ni r 
everywhere, we solve L k ' 1 Ct = r in part of the domain, 
where r = L k ‘’l£'V. This is closer to the continuous 
problem Lu = f. Here L is the continuous counterpart of 
the discrete operator L k l and u is the exart solution for the 
continuous problem; u k_1 is the exact solution to the discrete 
problem and 0 is the improved solution due to the 
modified right hand side of the discrete approximation. 

The deferred-corrector in (1) is input as a prescribed 
fimrtirtnal form on the right hand side of the fine grid 
equation. On any given grid level, the calculation is 
initially considered without the DC term. After a 
reasonable level of convergence is achieved, e.g., 10" , the 
DC term is evaluated from this known solution. This value 
is prescribed on the finest grid and added explicitly to the 
right hand side of the equations. This term is transferred 
to the coarser grid levels through the standard multigrid 
procedure. If the DC term is introduced earlier, divergence 
results. This is due to the fart that the solution is initially 
quite poor and therefore the prescribed DC term is 
significantly in error. This distorts the differential equation 
and induces an instability in the pressure during the 
relaxation procedure. If the RNS solution is allowed to 


converge moderately before introducing the DC term, the 
overall solution procedure, with the DC addition, exhibits 
no significant degradation in rate of convergence for the 
examples considered herein. 

For turbulent flow modelling, the eddy viscosity v, is 
calculated only on the finest grid. The fine grid values of 
v t are injected to the corresponding coarse grid points 
during the restriction step. The v, values, at points on the 
coarse grid, that lie outside of the extent of the fine gnd 
are not updated. This procedure is validated by 
calculations that do not include any adaptivity. In this case, 
the fine grid extent is the same as that of the coarse grid. 
Therefore, the v, values at all points on the coarse grid will 
be updated. The results obtained from this full refinement 
calculation and those from the fully adaptive multigrid 
calculations are identical. 

Tnterdomain T r ansfer of Boundary Cond i tions an d 
Conservation at Crid Interfaces 

For a given subdomain, the following boundary 
conditions are to be prescribed: 

u = v = 0 at q = 0; u = 1, p = 0 at r\ — rim*i> 
p t = 0 or p prescribed at f = 5^; u and v are given 
by freestream values at £ * £<>. 

For external flows, if a subdomain has its outflow at some 
l < then the boundary condition on pressure changes 
from Neumann to Dirichlet type. For internal flow, the 
outflow boundary condition is Dirichlet type for the 
pressure. Also, if the lower boundary of a subdomain is at 
some n > 0, then non-zero velocities have to be prescribed. 
In time-dependent, characteristic-based, Navier-Stokes 
computations, that use locally embedded grids, boundary 
conditions are required for all variables, i.e., u, v, and p. Tn 
addition, special care has to be taken to ensure that mass 
conservation is not violated locally or globally. 

In the pressure-based trapezoidal or box' formulation, 
this difficulty does not occur as the normal velocity v in n, 
or u in £, is not prescribed at the upper or lower, or 
outflow boundaries. Only the tangential component u is 
prescribed at the upper interface or interdomain boundary. 
The pressure-based box-type differencing allows for the 
calculation of the normal velocity at the outer boundaries 
and the pressure at the body surfaces. The normal velocity 
is computed from the continuity equation and therefore 
mass conservation is automatically satisfied on all levels, for 
all subdomains. This eliminates the need for special 
interpolation formulae to ensure conservation of mass when 
the boundary conditions are prescribed from the coarse grid 
solution. Thus weak instabilities, that arise when such 
methods are applied in Navier-Stokes formulations without 
satisfying conservation, do not appear in the present 
method. Direct evaluation of the pressure at the inflow or 
lower boundaries with the trapezoidal or box discretization 
also eliminates the need for special pressure boundary 
conditions. 

The calculation is performed sequentially rather than 
in parallel in the various subdomains. As such the 
boundary conditions at the inflow and outflow stations for 
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each subdomain are updated with the latest available 
values. The overlap allowed in the subdomaining process 
follows the following rules. 

a. The last station of any subdomain, which is at some 
l < Imu coincides with the first station of the 
subdomain to its right, (if one exists), where the 
pressure is computed. 

b. Similarly, the inflow station of any subdomain, which 
is at some £ > 0 coincides with the last station on the 
subdomain to its left, (if one exists), where the 
velocities are computed, 

c. If the inflow station or the outflow station of a given 
subdomain coincides with the physical boundaries of 
the global flow field then the boundary conditions 
discussed previously for inflow, outflow, upper and 
lower boundaries are used for these subdomains, 

d. If there are no subdomains to the right, for the cases 
in a), or if there are no subdomains to the left, for 
the cases in b), then these boundaries are updated 
using coarse grid values during the multignd 
prolongation process. 

In the vertical direction an implicit solver is applied and no 
overlap is necessary. 

If a subdomain has only one of its horizontal 
boundaries in common with that of another subdomain, 
then updating the boundary conditions along this edge, 
after one sweep in all subdomains, leads to iterative 
divergence on this subdomain. This influence gradually 
filters through to other subdomains. If these boundaries 
are updated through the multigrid transfer processes, then 
the calculation is convergent. This reflects the fart that an 
update of just one boundary after each sweep, with the 
other three updated only during the multigrid transfer 
process, leads to an inconsistency. This constrains the 
variables from adjusting to changes that occur dynamically, 
as the solution evolves in the various subdomains. The 
multigrid transfers provide the correct dynamic response to 
changes between subdomains. 

Results 

All of the calculations presented herein are initiated 
on the coarsest grid, with uniform flow velocity and 
pressure. On the finer grids, the interpolated coarse grid 
solution fields are, sequentially applied as initial 
approximations. Since convergence to the final solution is 
improved with more accurate initial approximations 5,9,10 , the 
adaptive multigrid framework introduces this element in a 
natural and convenient fashion. It is significant that for all 
of the examples presented herein, Reynolds number 
continuation is not required in order to obtain a solution 
for any of the prescribed values of Re. Even for highly 
interactive, large Re flows, the solution is obtained directly 
with uniform initial values at the designated value of Re. 

Example 1: Flow over a finite flat plate: the trailing 
edge problem 6 . 

Figure 1 depicts the grid obtained for a semi- 
coarsening (in £) adaptive calculation. The q grid is highly 
stretched and fixed. Note that the finer grids zoom in 


around the trailing edge, which is located at £ * 1.0 (the 
figure is scaled by a factor of 2 in the £-direction). 
Significantly, the extent of the finer grids in the q direction 
is progressively reduced, even when adaptivity is specified 
only in the £ direction. Although each multigrid level 
contains only one subdomain (in £) that requires further 
refinement on subsequent levels, the q extent of this 
subdomain is affected. 

Figure 2 depicts the composite grid obtained with full 
two-dimensional adaptivity. Within each subdomain, 
uniform grids, in both the £ and q directions, are 
prescribed. Figure 2 is an overlay of seven multigrid levels, 
each of which comprises several subdomains. In each level, 
it is found that the subdomain, for which refinement in 
both directions is required, is centered around the trailing 
edge. The adaptive computations, both semi-coarsening 
and two-dimensional, are compared with non-adaptive 
semi-coarsening multigrid calculations. For the latter, a 
uniform fine grid in £ and a highly stretched q grid is 
prescribed. The grid stretch factor is chosen from the 
specified minim um and maximum Aq values, and the 
location of q^ that was applied for the two-dimensional 
adaptive study. The same q grid is employed for the 
adaptive semi-coarsening calculation. Figure 3 shows a 
comparison of the pressure coefficient Cp for the three 
calculations. There is good agreement in the pressure 
variation and, in particular, the predicted peak pressures. 
Table 1 summarizes the computer memory and CPU 
requirements. These are given as percentages of the non- 
adaptive, semi-coarsening, calculation. Note that the 
memory requirement for the one- and two-dimensional 
adaptive calculations are similar. This signifies that the 
specified q stretching for the semi-coarsening calculation is 
reasonable. 


Table 1. Summary of Computer Resource Requirements 
for the finite flat plate calculation 



Two-D 

Adaptive 

One-D 

Adaptive 

Full Refinement 

Aspect 

with stretched 
h grid 

CPU 

18.03% 

15.10% 

100.0% 

Memory 

12.90% 

13.22% 

100.0% 


The adaptive grid of Figure 2 defines the extent of 
the interaction zone surrounding the trailing edge. From 
large Re asymptotic triple deck theory, three layers with 
different length scales have been identified 12 , viz., a lower 
viscous rotational deck of 0(Re’ 5 ^), a middle inviscid 
rotational deck of (KRe" 4 ^), and an upper inviscid 
irrotational deck of 0(Re' 3/8 ). Since the vorticity is zero in 
the upper deck, and since vorticity is the monitored 
parameter for refinement in the q direction, the adaptive 
procedure should lead to a grid that does not require q 
refinement in this 'upper deck' region. The grid obtained 
from the two-dimensional adaptive calculation displays this 
result quite clearly. At each multigrid level, there is a 
region away from the body that is in fact refined only in the 
£ direction. This region, in the finest multigrid level, 
represents the extent of the upper inviscid irrotational 



region. Estima tes for the extent of the other two decks 
are also obtained from the grid structure. In more 
complicated flows, e.g., turbulent flow past the same 
geometry, for which analytical methods cannot be easily 
developed, the appropriate resolution in each distinct 
region will be automatically captured with the present 
multigrid adaptive procedure. In this sense, the 
computation results in a form of discrete asymptotic analysis. 

Example 2. Flow over a trough 6 . 

The second geometry to be considered is the laminar 
flow over a trough configuration. Both unseparated and 
reverse flows are computed with the two refinement 
strategies previously discussed. The trough surface is 
specified by y b = -D sech[4(x-Xo)], where D represents the 
maximum depth, which occurs at the location Xq. The 
values x<,=23 and Re - 80000 are used for the present 
calculation. The grid obtained from the two-dimensional 
adaptive procedure, for D=0.03 and with a region of flow 
reversal, is shown in Figure 4. Note that refinement m the 
n direction extends to a significantly greater distance than 
was found for the trailing edge geometry. This is due to 
the fact that the maximum vortidty now occurs near the 
outer edge of the separation bubble and not at the surface. 
Also note the sudden increase in the extent of the region 
where q refinement is performed. This signifies the 
increase in boundary layer thickness as a result of flow 
separation. The reversed flow region is essentially vorticity 
free; however, the current refinement strategy assumes that 
regions requiring refinement, in the q direction, will always 
have a lower boundary at the wall. This condition can be 
modified to allow for multiple q subdomains in the 
recirculation region. This was not considered necessary for 
the current calculations. Figure 5 depicts the pressure 
variations obtained for the three calculations discussed 
previously for purely attached flow and D= 0.015. Figure 
6 provides comparisons of the skin friction for the 
separated (D=0.03) case. Once again good agreement is 
obtained and significant gains in computer resource 
requirements are found (Table 2). The locations of the 

Table 2. Summary of Computer Resource Requirements 

for the trough geometry 


Aspect 

Geometry 

Two-D 

Adaptive 

One-D 

Adaptive 

Full 

MG 

with non 
uniform 
q grid 


Trough 

18.03% 

— 

100.0% 

CPU 

i 

(Sep) 

Trough 

(U0«p) 

7.10% 

16.80% 

100.0% 


Trough 

1632% 

— 

100.0% 

Memory 

(Sep) 

Trough 

(Uneep) 

5.10% 

63.4% 

100.0% 


separation and reattachment points computed by the two- 
dimensional adaptive calculation are at f =231 and 
l *254, respectively; the values predicted by non-adaptive 
full multigrid refinement are at 5 =231 and 5 =233. This 


further confirms the validity of the domain decomposition 
approach and the advantages of adaptive multigrid over full 
multigrid. All results presented here are in excellent 
agreement with all earlier results 1 -’ presented for these 
geometries. 

Example 3. Internal flow in a back step channel: 
laminar 73 and turbulent flows. 

For this flow, which is dominated by rather large 
recirculation regions, it is still possible to cany out the 
calculation for all Re considered herein by prescribing 
uniform initial flow conditions. This is true even for the 
relatively difficult, although somewhat artificial two- 
dimensional calculation with Re =800 (based on channel 
height). For this Re value, two separation bubbles, one on 
each wall are evident. Reynolds number continuation, as 
applied in many other reported NS solvers , is still not 
required for the present calculations. Both laminar 
(Re =800) and turbulent results are obtained for the back 
step geometry. The standard two equation k-e model 
discussed earlier is used for turbulence closure. 

For laminar calculation, a step height to charnel 
height ratio of 03 is used. The reattachment length (Xr) 
for the primary recirculation zone is compared in Table 3 
for a range of laminar Reynolds numbers. Comparisons 
are given for the present 2-D adaptive method, full 
refinement with the standard non-adaptive mulhgnd 
method, and earlier calculations by Ferager 5 , Caruso and 
Sodropoulos 13 . The calculated reattachment length for the 
adaptive and non-adaptive procedures is identical to two 
decimal places; however, the computational effort is 
considerably less for the former, see Table 4. 

Table 3* Comparison of Reattachment Length for 
Laminar Backstep Channel Flow 


j 

Present Calculations 

PPP 

REF 

[15] 

REF 

[13] 

I Re 

Adaptive 

Non- 

Adaptive 

xvnr 

[5] 

1 133 

1.94 

1.94 

2.0 

1.95 

1.84 


335 

335 


335 



432 

432 


4.40 



530 

530 


5.40 



For adaptive refinement in the q direction the 
truncation error is scanned from the wall towards the outer 
boundary. For external flow, the vorticity gradient 
decreases exponentially and a thin layer near the wall, 
where refinement is maximum, can be identified. This 
region is specified by fixing the upper boundary at the 
furthest point, or largest r\ value taken over all 5 locations, 
that satisfies the truncation error tolerance. For uternal 
flows, boundary layers, where refinement in q should be 
required, exist at both boundaries in the normal or q 
direction. However, the number of grid points that are 
necessary to resolve the flow gradients in the q direction is 



quite moderate and therefore no attempt was made to 
adaptively refine in this direction. Instead, the full 
multigrid procedure is applied in the q -direction for each 
of the subdomains for which ^-refinement is necessary. 
This allows for different uniform q grids in the different l 
subdomains. 


Table 4. Summary of Computer Resource Requirements 
for the Backstep Geometry 


Aspect 

Re 

Adaptive Multigrid/Full 
Refinement NonAdaptive Multigrid 


133 

35.49 % 

CPU 

267 

36.15 % 

400 

4623 % 


600 

50.40 % 


133 

30.80 % 

Memory 

267 

37.44 % 

400 

41.49 % 


600 

47.49 % 


Table 4 displays the computer resource requirements for 
the backstep channel calculation. For each Reynolds 
number, the CPU and memory requirements are shown as 
percentages of the corresponding non-adaptive calculations. 
Note that as the Reynolds number increases from Re * 133 
to Re =600, the number of grid points required to resolve 
the flow field increases. This is expected, as the size of the 
separation bubble increases with Reynolds number. The 
number of required grid levels, as well as the finest mesh 
size for all Reynolds numbers up to Re =600, is identical in 
each direction. A total of five multigrid levels are defined 
for all Reynolds numbers up to Re =600. However the 
subdomain extent for each multigrid level is different for 
different Re. The extent of the finer grids is governed by 
the size of the recirculation zone, which increases as Re is 
increased. For the Re =800 case, six multigrid levels are 
required, as the change in the solution from level 4 to level 
5 is still significant and greater than the specified tolerance. 

An increase in computational time and memory 
requirements is observed as the Reynolds number and 
associated number of grid points increases. The time 
required for the full refinement non-adaptive calculation 
increases only marginally as Re is increased from 133 to 
600. This increase is entirely due to the changing nature of 
the flowfield, since the same grid is used throughout this 
Re range. More specifically, when the degree of velocity 
relaxation is increased due to the increasing extent of the 
reversed flow region, the convergence rate degrades and 
additional iterations are required to achieve the specified 
tolerances. Furthermore, the percentage gain in adaptive 
over non-adaptive procedures is reduced as the Reynolds 
number increases, e.g., to about 50% at Re =600. 

The effect of location of the outflow boundary and the 
non-reflectivity of the outflow boundary conditions are 
important aspects of this study. The adaptive multigrid, 
do main decomposition procedure is initiated on a very 
coarse grid, and yet, it is possible to place the outflow 
boundary quite far downstream, e.g^ 60 step heights, and 
still recover very accurate and computer efficient 


computations. The solutions at the outflow are in almost 
perfect agreement with the analytic fully developed flow 
values. Although the finer grids in l and n occur in 
subdomains much further upstream, in and near the reverse 
flow regions, the influence of the outflow boundary 
conditions is propagated through the multigrid transfers to 
and from the coarser grids connecting the various 
subdomains. This allows for efficient transfer of 
information without excessive grid specification. In 
addition, as is shown 8 , the RNS pressure flux-splitting 
procedure allows the outer boundary to be placed very far 
upstream, e.g., within the upper wall recirculation region, 
without solution degradation. 

For the laminar backstep, the DC terms neglected in 
the RNS approximation have been included after obtaining 
a reasonably converged base solution for the RNS system. 
For this geometry, the vertical wall region near the step 
corner represents the only portion of the flowfield where 
the full Navier Stokes terms are of any consequence. 
Along the vertical wall, the v„ term represents the vortical 
or diffusive boundary layer influence. It is found 8 that in 
this region, although the inclusion of the DC term does not 
produce a significant quantitative difference, some 
qualitative difference is observed in the solution. Figures 
7a-7d depicts the streamwise velocity profile for Re =400 at 
four successive stations near the corner. Note that the 
effect of the DC diminishes rapidly away from the corner. 
A significant difference between the two solutions is 
associated with a small positive axial velocity that is 
observed near the corner when the DC is included. This 
represents a counter rotating vortex within the primary 
separation bubble. The reattachment length remains 
unchanged even when the DC is included. 

The non-reflectivity of the Dirichlet pressure outflow 
boundary condition was tested for the severe Re =800 case 
with calculations on computational domains of three 
different lengths. It was found that for the cases 
considered, locating the outflow boundary further upstream 
did not have a significant effect on any of the solutions. 
Comparisons of the streamwise velocity profile at the 
streamwise location X=7 are given in Figure 8. The 
outflow boundary was located at X = 7, X=15 and X=30; 
the results clearly indicate that the effect of the outflow 
location on the solutions is minimal Figure 9 shows the 
comparison of the stream function contours obtained by 
placing the outflow at X=7 and X=15. The two contour 
patterns are identical 8 . 

A benchmark solution for Re =800 has been 
published by Gartling 14 . The present solution, which is 
obtained with the adaptive multigrid domain decomposition 
procedure, is compared with these results. In Figure 8, 
comparison with the benchmark solution of the streamwise 
velocity profiles at X*7 is also shown. Note that reverse 
flow also occurs on the upperwall. The appearance of this 
upper separation bubble is thought to introduce three- 
dimensionality into the flow, and for this reason there is 
some disagreement between the experimental results and 
all of the numerical solutions. However, the present 
results, which are totally grid independent, agree quite well 
with most of the other numerical computations. Due to the 
very fine meshes that have been prescribed with the 



multigrid domain decomposition procedure, the residuals 
and truncation errors are quite small and therefore these 
numerical solutions, are considered to be highly accurate. 
Figure 10 depicts comparison of vorticity profiles. The 
agreement with the benchmark contours 14 is excellent 
throughout Moreover the results are essentially unchanged 
with the outflow boundary at X=30 or X=7, which is inside 
the upper recirculation zone 8 for the Re = 800 case. 

Figures 11 and 12 show typical grids obtained using 
the adaptive multigrid procedure for the backstep channel 
Figure 11 shows the grid for Re = 100 and Figure 12 depicts 
the grid for Re =267. Note that the region covered by the 
finest grid is larger in extent for Re =267. This is expected 
as the region of reversed flow is larger in this case. Note 
that the grid used in the y-direction is quite coarse far 
downstream where the flow is fully developed. This reflects 
the fact that the truncation error is very low in the normal 
derivatives in this region. Since the fully developed flow 
profile is parabolic, central differencing will theoretically 
incur zero truncation error. In streamwise direction, a grid 
as coarse as Ax =0.5 can be used towards the outflow 
without loss of accuracy. The adaptive multigrid procedure 
clearly utilizes this fact and provides optimal resolution. 

The turbulent flow past the backstep channel is of 
interest in many engineering applications. For the current 
study, the k-« model has been used to compute the 
Reynolds stress terms. This model requires that inflow 
values for k and e be prescribed. For the present 
calculation these values are generated from a straight 
channel turbulent flow computation with a step initial 
profile. The k, e, u and v profiles obtained at the outflow 
of the straight duct are used as the inflow conditions for 
the backstep calculation. This inflow is located at a 
distance five step heights upstream of the step comer. It is 
noted that the overall nature of the flow field, including the 
reattachment length is strongly dependent on the inflow 
values used for k and e. Different profiles for k and e can 
be generated by varying the length of the straight duct. 
This will greatly alter the backstep channel solution. A 
channel length of 23 step height was used to generate the 
inflow profile for the present calculation. This ensures that 
the profiles are quite well developed. The corresponding 
velocities, k an d e profiles are then used for the backstep 
channel calculation. Once again, the inflow is located five 
step heights upstream of the step comer. Thangam et aL 
observed that the predicted reattachment length using the 
standard isotropic k-e model was in error by about 12%, 
when compared with the experimental value. They showed 
that a modified k-e model which takes anisotropic effects 
into account provides improved results. In another study by 
Thangam 17 , it was shown that a modified isotropic k-e 
model can also lead to improved results. This model, 
which requires only the variation of a single constant, is 
considered for some of the calculations presented herein. 

The step height to channel height ratio used is 1-3 in, 
the present calculations. A reattachment length (x/H) of 
7.04 was obtained for a Re =132000 based on channel 
freight. This is in very good agreement with the 
experimental value of 7.1. Figure 13 shows the 
streamfunction contours for the same flow. Note that there 
is a secondary counter rotating vortex within the primary 


separation bubble. This was also observed by Thangam et. 
al . The calculation was also performed by modifying the 
value of c «2 to 11/6 as suggested 17 . The reattachment 
length increased, in this case, to 7.66. This was also the 
trend observed 17 . Although in their case, the standard k-« 
model was underpredicting the reattachment length and 
modifying the constant produced acceptable results. 

Summary 

An adaptive multigrid domain decomposition method 
has been used to efficiently compute incompressible 
laminar and turbulent flows with regions of recirculation 
and strong pressure interaction. A low order RNS system 
of equations, a fully consistent primitive variable non- 
staggered grid solver, accurate mass conservation at 
subdomain interfaces and global boundaries, non-reflective 
outflow boundary conditions and a pressure-based flux-split 
discretization are the key features of the procedure. The 
adaptive multigrid domain decomposition procedures allows 
for efficient grid definition consistent with asymptotic 
theory and for effective transfer of information to and from 
fine grid high gradient regions to coarse grid 'inviscid' 
regions. 

significant gains in computer resources have been 
achieved when compared to standard non-adaptive 
methods. Good agreement is obtained between the present 
solutions, standard non-adaptive full refinement 
computations and other published results. The 
computational cost is several times smaller than that 
required by most other NS methods 14 . For example, the 
CPU required for the backstep c han nel calculations, with 
Reynolds numbers in the range 100-600, varies between 5- 
10 minutes on an IBM 320 RISC/6000 workstation. For 
the Re =800 case, an additional multigrid level is added to 
ensure grid independence and the CPU required is 
increased to approximately 30 minutes on the same 
workstation. All solutions are initiated with uniform flow 
approximations and Reynolds number continuation is not 
required, even for the relatively complex Re =800 case. 
Grid convergence has been established efficiently through 
an adaptive multigrid procedure. The outflow boundary 
condition has been shown to be non-reflective. In addition, 
it has been shown that the procedure is not very sensitive 
to the location of the outflow, i.e. far downstream or 
somewhat closer to the inflow. The flux-split discretization 
allows for direct computation of the normal velocity and 
therefore mass conservation at grid interfaces and 
subdomain boundaries is achieved in a simple fashion. 
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